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Abstract 


This is the forth report in a series of technical reports that describe 
separated two-phase flow model application to the cryogenic loading op- 
eration. In this report we present the structure of the code. The code 
consists of five major modules: (i) geometry module; (ii) solver; (iii) 
material properties; (iv) correlations; and finally (v) stability control 
module. The two key modules - solver and correlations - are further 
divided into a number of submodules. Most of the physics and knowl- 
edge databases related to the properties of cryogenic two-phase flow are 
included into the cryogenic correlations module. The functional form of 
those correlations is not well established and is a subject of extensive re- 
search. Multiple parametric forms for various correlations are currently 
available. Some of them are included into correlations module as will be 
described in details in a separate technical report. Here we describe the 
overall structure of the code and focus on the details of the solver and 
stability control modules. 


1 Introduction 

The algorithm of the separated two-phase flow (see also [LuchDG-I]) 
originates from Liles and Reed [Liles-78], which is in turn based on Har- 
low and Asden [Harlow-68] all-speed “Implicit Continuous-Fluid Eule- 
rian (ICE)” algorithm. In the current work we implemented the nearly- 
implicit extension of this algorithm following closely the results of [RELAP5- 
I]. The ICE-based algorithms is usually called “weakly-compressible”, 
even though the compressibility effects are fully accounted for in numer- 
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Eigure 1. A simplified structure of the code, (top) Eour main blocks of 
the algorithm. Bottom left: Structure of the boundary conditions block. 
Bottom right: Eour main sub-steps of the integration step. 
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ical discretization. Almost all reactor thermal-hydraulics codes belong 
to this class and utilize first-order finite-difference donor-cell/upwinding 
based schemes, implemented on structured staggered meshes, see e.g. [Ran- 
som89, TRACE, Nourgaliev] for further discussion In what follows we 
present the details of the algorithm as it was codded. A simplified struc- 
ture of the code is shown in the Fig. 1. The algorithm consists of the 
four main blocks 

1. geometry, 

2. initial conditions, 

3. boundary conditions, 

4. one integration step. 

The geometry and initial condition blocks are set once in the beginning 
of the code execution. The last two blocks are repeated in the loop N- 
steps until integration is completed. These two blocks have a complex 
structure. 

The boundary conditions block includes calculation of the interphase 
geometry, the heat and mass fluxes at the liquid/gas/wall interphases, 
mass and energy fluxes through the dump valves and through the in- 
put/output valves. This block includes the pressure drop and heat trans- 
fer correlations for the two-phase flow. These correlations are based an 
the flow patterns recognition and have very nontrivial structure that will 
be described in details in a separate technical report. The integration 
block consists, in turn, of four main sub-steps: 

1. first step of integration, 

2. second step of integration, 

3. control of the values of dynamical variables, 

4. time step control. 

The first two sub-steps are designed in the spirit of the predictor- 
corrector architecture of the nearly-implicit method. The last two sub- 
steps are developed to deal with multiple instabilities inherent to the two- 
phase flow algorithms and are essential for the smooth code execution, 
see e.g. [Nourgaliev, Cordierl3, LuchDG-I]. 

Out of four integration sub-steps the first one has the most nontriv- 
ial structure and its details are the key to the successful performance of 
the algorithm as a whole. It involves calculations of the upwind, face- 
centered, and volume-centered values of the dynamical variables. The 
overall stability and accuracy of the algorithm in the quasi-equilibrium 
limit (no heat transfer) is mainly determined by the solution of the two 
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matrix equations: (i) matrix equation for the expanded form of the con- 
servation laws for the phasic masses and energies and (ii) matrix equation 
for the phasic velocities (see Fig. 1). 

Before we discuss the structure of each block in more details, we 
remind for convenience the model equations of the nearly-implicit scheme 
method. 


2 Model 


2.1 The model equations 

We model cryogenic loading and chilldown using Wallis equations [Wal- 
lis69] for a one-dimensional separated two-phase flow. The model con- 
sists of a set of conservation laws for the mass, momentum, and energy 
for the gas 

{AaPg)^ + {AaPgUg)^^ = ATg 
{AapgUg) ^ -h (^Aapgufj ^ + Aap^x = -AapgZ^x 
Tgw^wg Tgi^i T -^^g^ig 

{AapgEg) ^ + (AapgEgUg)^^ = Attp - (pAaUg)^^ 

T Qgw^wg T Qgi^i T -^^igHig AE yjgHyjg 

and for the liquid 


{Aj3piui)^ + {A(3piuf)^^ + APp^x = -A(3piz^x- 

'^Iw^wl AE gll^l (2) 

{ApEipi) ^ + {Al3Eipiui)^^ = A/3p^t ~ {pAj3ui)^^+ 

Qlw^wl T Qli^i igHig AE yjgHyjg 


coupled to the equation for the wall temperature 
dT 

PwC-wdw — h^g (Eg lui) -|- h^i (El Tiu) -|- hamb {Eamb l~u)) . (3) 


Here p, a, T, and p are pressure, temperature, and density of the fluid. 
E is the total specific energy, Hig and H^g is the specihc enthalpy of 
the gas generated at the interface and near the wall respectively, u is 
the fluid velocity and h is the heat transfer coefficient. Other model 
notations are explained in the Section Nomenclature. 


2.2 Discretization 

The equations (1), (2), and (3) are integrated on a one dimensional grid 
shown in the Fig. 2. The energy and density conservation equations are 
integrated over N control volumes centered at locations L = 1 , . . . , 
shaded by yellow. The momentum equations are integrated over N — 1 
control volumes of the staggered grid centered at locations j = 1, ... ,N— 
1 shaded by blue. The equations for the wall temperature are integrated 
over N control volumes shaded by green. 
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Figure 2. Grid grid with N control volumes shaded by yellow, staggered 
grid with N — 1 control volumes shaded by blue, and N control volumes 
for the wall shaded by green. 

2.3 Brief summary of the nearly implicit scheme 

The algorithm described below is a variation of the nearly-implicit method [RELAP5- 
I] and [RELAP5-VI] (see also [LuchDG-I]). 

2.3.1 Boundary conditions 

Boundary conditions module includes analysis of the interface geometry, 
phasic masses, heat fluxes, and flow through the valves. At present 
this module involves a simplified set of pressure drop and heat transfer 
correlations based on the effective geometry of conceptually stratified 
flow (i.e. liquid and gas are assumed to be always stratified). The 
calculations within this module can be briefly summarized as follows 

• determine geometry of the phasic interface; 

• calculate frictional losses; 

• find heat transfer coefficients at the wall; 

• find mass and enthalpy fluxes at the liquid/gas interface; 

• calculate mass and enthalpy fluxes through the dump and input / output 
valves. 

The extended version of the boundary conditions module will include 
pressure drop and heat transfer correlations based on the flow pattern 
recognition and will be discussed in details in a separate technical report. 

2.3.2 First sub-step 

The calculations of the velocities, pressure, and provisional values of the 
masses and energies at the first sub-step of the algorithm is the key to 
the stable performance of the nearly-implicit method. These calculations 
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are designed to break limitations of the material CFL and to increase 
implicitness of the method. This sub-step is structured as a predictor 
of the fractional time step technique. The CFL limitations are lifted by 
implicit calculations of the new velocities and pressures in the system. 
The implicitness is further increased by estimations of the provisional 
values of the densities, energies, heat and mass transfer coefficients. The 
calculations can be briefly summarized as follows; 

• Calculate upwind, face-centered, and volume-centered velocities; 

• Solve expanded equation with respect to pressure in terms of new 
velocities; 

• Substitute this solution into momenta equations and solve resulting 
block tri-diagonal matrix equation for the new velocities; 

• Find new pressure; 

• Find provisional values for energies and void fractions using ex- 
panded equations; 

• Find provisional values of mass fluxes and heat transfer coefficients 
using provisional values of temperatures obtained. 

Corrections to the provisional values of void fractions, densities, and 
energies are found at the second sub-step of the integration step of the 
algorithm. 

2.3.3 Second sub-step 

To hnd new (corrected) values values of the densities, void fractions, 
and energies we solve unexpanded conservation equations for the phasic 
masses and energies using fully implicit method. The solution is reduced 
to independent solution of four tridiagonal matrices. The values of pres- 
sure and velocities in these matrices are taken at the new time step. 

2.3.4 Control 

Despite the enhanced stability of the nearly-implicit method the insta- 
bilities remain a challenging problem. Multiple sources of instabilities 
in modeling separated two-phase flow include non-hyperbolicity, lack of 
positivity, and phase appearance/disappearance (see e.g. [LuchDG-I] and 
references therein). To obtain a reliable and stable execution of the code 
extended controls of the dynamical variables and of the time step are 
introduced into the code. These controls include 

1. check that void fraction is between 0 and 1; 

2. if void fraction is within the predefined range of the boundary 
values corresponding to the phase appearance/disappearance apply 
smoothers to the temperatures, densities, and velocities; 
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Figure 3. A set of control volumes introduced for the transfer line model 
of the KSC testbed. 

3. check if phasic temperatures are within thermodynamic range; 

4. check pressure range; 

5. check mass conservation in each control volume; 

6. check mass conservation in the system as whole; 

7. compare the provisional and hnal values of phasic densities and 
energies; 

8. if the dynamical variables are outside of the physical range or vari- 
ation if variables is too large the integration step is repeated with 
reduced value of the time step At = At/2, while At > At^m- 

Once control checks are completed the calculations return to the 
step 2.3.1. The sub-steps 2.3.1 through 2.3.4 are repeated until final 
time of the integration T = is reached. 

3 Geometry Module 

One of the important properties of the method briefly outlined above is 
the ability to resolve the two-phase flow dynamics and heat transfer for 
a complex geometry of the pipes. The geometry of the model is defined 
as an ordered collection of N control volumes Ai = {CV}^ . The finite 
volume mesh Ai is characterized by the following constant parameters: 
(i) length Axi, (ii) perimeter (iii) height hi, and (iv) inclination angle 
Oi. The following related geometrical parameters are also frequently 
used: (v) diameter di, (vi) surface area Si, (vii) volume Vi, and (viii) 
cross-sectional area Ai. 
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Figure 4. A single control volume of the model. 


For example, the control volumes introduced for the cryogenic testbed 
are shown in the Fig. 3. The volumes are numerated at the bottom of 
the figure. The blue shaded areas indicate the quantity of liquid in each 
control volume. A single control volume of the model is shown in the 
Fig. 4. A number of important geometrical parameters of the control 
volume that are shown in the figure change every time step and belong 
to the set of dynamical variables of the system. These parameters are: 
(i) dry perimeter lgw,i, (h) wetted perimeter kw,i, (hi) interface perimeter 
li^i, (iv) height of the stratified liquid level bi^i. 

3.1 Components 

In addition to the standard geometrical parameters there is a large num- 
ber of components that dramatically affect the fluid flow and have to 
be taken into account in any industrial cryogenic loading system. These 
components include e.g. valves, bends, hlters, httings, pumps etc. The 
components are characterized by their location and frictional and heat 
losses. For example the locations of the valves that control flow in cryo- 
genic testbed at KSC together with location of pressure and temperature 
sensors are shown in the Fig. 5. The solid black line corresponds to the 
flow path. The location of the components including dump valves, con- 
trol valves, and sensors can also be resolved in the figure. 

These information is integrated into the code in the form of tables. 
An example of the table is shown in the Fig. 6. 

In this way the information about components can be attributed 
to each control volume for arbitrary set of control volumes A4. For 
example the location of the dump valves in the set of control volumes 
corresponding to the KSC testbed model is shown at the top of the Fig. 3. 

The details of the calculations of the source terms that include in- 
formation about pressure and heat losses in the system components are 
given the following section. 
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4 Boundary Conditions Module 


Boundary conditions within this report are interpreted in a more general 
sense than usual. They include not only conditions at the input/output 
valves, but also conditions at the liquid/gas/ wall phasic boundaries and 
at a liquid/gas interface. Therefore, the current version of the boundary 
condition module includes calculations of the simplified correlations rela- 
tions. However, the future versions of the algorithm will have a separate 
correlation module. 

Note, that the boundary conditions for the continuity and energy 
equations are determined on the grid of the control volumes of the sys- 
tem, while boundary conditions for the momenta equations are dehned 
on the staggered grid. To this end it is convenient to rewrite equations 
(1), (2) in a matrix form 


dU OF _ dV dG_ 

dt dx dt dx 

Here the conservative variables are 


U 


Aapg 

AapgEg 

APPI 

A/3piEi 


AaPgUg 

Al3piui 


( 4 ) 


( 5 ) 



Figure 5. Geometry of the flow path is shown by the solid black line. 
The dump valves (DCV) are indicated by yellow circles. The control 
valves (CV) are shown by cyan circles. The temperature sensors (TT) 
are shown by bike circles and the pressure sensors (PT) are shown by 
red circles. 
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Figure 6. An example of the table with a list of cryo-testbed components 
integrated into the model. 


and fluxes are 

AapgUg 

P ^ AaUg {PgEg + p) 

AI3piui 

AI3ui {Eipi + p) 

The Boundary Conditions Module is primarily concerned with the 
calculation of the source terms 


G = 


AapgV? 

Al3piuf 


( 6 ) 


Su 


ATg- rh^i 

Aotp^t T Qgwlwg T Qgi^i T -^^igEig + AT'^gH^g TTl^ihg 

-AFg 

T Qlw^wl T Qli^i igHig AFyjgliy^g 


( 7 ) 


Eo,p^x AopgZ^x T~gw^wg Tgi^i T EYgUig 

^ +Al3p^x - Al3piz^x - Tiiulwi - Tiik - ATgUil 

The source terms (7) are calculated twice during one integration step. 
Once during the Hrst sub-step and once during the second sub-step. 
The source terms (8) are calculated only once before the solution of the 
velocity matrix equation. 

The source terms are calculated using a few sub-steps. The first 
sub-step involves the following calculations: 

• determine flow patterns for each control volume; 

• calculate geometrical parameters of the flow based on the type of 
the flow pattern; 
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• find heat fluxes at the liquid/gas/wall interfaces for the values of 
thermodynamic variables found at the previous time step; 

• determine the mass and enthalpy fluxes using values of the heat 
fluxes and latent heat of vaporization/condensation; 

• find mass/enthalpy fluxes through the dump valves and input /output 
valves. 

Once the terms Sy are calculated, the new values of pressure are 
found by solving the expanded and modihed version of the first equation 
in (4) with respect to pressure as described in more details in Sec. 5.1. 

At the next sub-step the new values of pressure and old values of 
the thermodynamic variables are used to find source terms in (8). The 
values of Sy obtained at this sub-step are used to find the new velocities 
by solving the non-conservative version of the second equation in (4) as 
discussed in details in Sec. 5.2. The new velocities and new pressures 
are used to find provisional values of the thermodynamic variables (see 
Sec. 5.3). 

Next, the values of the Su are updated using provisional temperature 
and mass fluxes and the conservative form of the hrst equation in (4) is 
solved to find new values of thermodynamic variables as described in 
Sec. 6. 

It can be seen from this brief description that all the non-trivial 
physics of the problem and related calculations are coupled with an anal- 
ysis of the source terms. Further details of this analysis are discussed in 
the following few subsections of this section. 

4.1 Flow patterns 

Mass and heat transfer between the phases of the two-phase flow is closely 
related to the geometrical patterns formed by the flow. Examples of 
typical patterns include [BejanOS] 

Stratified flow. At low liquid and gas velocities, there is complete sep- 
aration of the two phases, with the gas in the top and the liquid in 
the bottom, separated by an undisturbed horizontal interface. 

Bubbly flow. In this regime, the gas is dispersed in the form of discrete 
bubbles in the continuous liquid phase. The shapes and sizes of 
the bubbles may vary widely, but they are notably smaller than 
the pipe diameter. 

Slug flow. Increasing the gas fraction, bubbles collide and coalesce to 
form larger bubbles similar in size to the pipe diameter. These have 
a characteristic hemispherical nose with a blunt tail end, similar to 
a bullet, and are referred to as Taylor bubbles. Successive bubbles 
are separated by a liquid slug, which may include smaller entrained 
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Figure 7. Example of the flow patterns in the two-phase flow. Nota- 
tions for the flow boundaries are (following [BejanOS]): (A) annular flow; 
(D) dispersed flow (droplets flowing in the gas); (I) intermittent flow 
(switches between patterns of slug and annular flow); (M) mist flow; (S) 
stratified flow; (SW) wavy stratified flow. 


bubbles. These bullet-shaped bubbles have a thin film of liquid 
between them and the channel walls, which may flow downward 
due to the force of gravity, even though the net flow of liquid is 
upward. 

Annular flow. Here the bulk of the liquid flows as a thin film on the 
wall with the gas as the continuous phase flowing up the center of 
the tube, forming a liquid annulus with a gas core whose interface 
is disturbed by both large-magnitude waves and chaotic ripples. 
Liquid may be entrained in the high-velocity gas core as small 
droplets; the liquid fraction entrained may be similar to that in 
the film. This flow regime is quite stable and is often desirable for 
system operation and pipe flow. 

Mist flow. When the flow rate is increased even further, the annular 
film becomes very thin, such that the shear of the gas core on 
the interface is able to entrain all the liquid as droplets in the 
continuous gas phase (i.e., the inverse of the bubbly flow regime). 
The wall is intermittently wetted locally by impinging droplets. 
The droplets in the mist may be too small to be seen without 
special lighting and/ormagnification. 

The detailed analysis of the flow patterns goes beyond the scope 
of the present discussion and will be given in a separate technical re- 
port [LuchDG-III]. Here we only provide an example of how flow pat- 
terns are determined in practice. Following the results of Kattan et 
al [Kattan:98a] the boundary for the stratified flow on the plain (m, x) 
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we have 


’^strat 


{226.3fAi^dAl^aPG {pL - Pg) f^L9 ] 

Y^n-ylvr^ ' 


1/3 


+ 20x, 


(9) 


for the wavy flow one can obtain 


m. 


^^^IwddPPLPG 


wavy 


7T 

256f 


- {2bi - 1)^] 
and for the mist flow the boundary is 




^mist — 


'iQSQA^^^^gDpLpG 

^Fr\' 


KWeJL_ 


.5 


( 10 ) 


( 11 ) 


Here the following parameters were introduced. The friction factor 


^Ph 


(l.l38 + 21og 


) 

\.5A[^d J 


-2 


The Weber number Wer, = — — , the Froude number Frr = , 

^ PLcr ’ ^ Pl 3^ 

the gas quality y = ■ "W , and the departure from nucleate boiling 

ilL'u 

Qdnb = OASIpq"^ H ig[g {pL — Pg) cr]^^'^ ■ The two fitting functions Fi 
and F 2 introduced in (9) - (11) have the form 


Fi (q) = 646.0 


q 

qONB 


2 

+ 64.8 


q 

qONB 


and 

-^ 2 ( 9 ) = '^^■^{q/qDNB) + 1.023. 

We can see that a number of the geometrical parameters of the flow 
has to be determined before the boundaries between the flow pattern 
can be calculated. These parameters include e.g. cross-sectional area of 
the gas flow normalized by the pipe diameter Ag^d, height of the liquid 
flow bi, etc.. We now describe how to calculate geometrical parameters 
of the flow. 

4.2 Geometrical parameters of the stratified fiow 

To find geometrical parameters of the flow one assumes usually [BejanOS] 
that the flow is “conceptually” stratified for all temperatures and flow 
rates, i.e. the flow cross-section has the form shown in the Fig. 8. 

The following parameters of the stratified flow have to be determined 
(see Fig. 8). First, the stratification angle 9 is found by noticing that 
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Figure 8. Geometry and parameters of the stratified flow in the pipe 
cross-section. 

the liquid cross-section area Ai = (1 — a) A (shaded by blue color in the 
flgure) is related to 9 as follows 


Once this equation is solved with respect to 9, all other required geomet- 
rical parameters, including 

• bi - height of the liquid level, 

• li - perimeter of the interface, 

• Igw - perimeter of the dry wall, 

• liui - perimeter of the wetted wall, 

• Ag - cross-section area of the gas, 

• A[ - cross-section area of the liquid, 

• Sg - dry area of the wall, 

• Si - wetted area of the wall, 

are found using simple geometrical relation. 

Once the equation for the 9 above is solved for a given void fraction 
a the effective liquid level height and other geometrical parameters can 
be found as follows: 


(1 — a)A = ~ sin0). 
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^iw — itD 


gw, 


li = D sin 



Similar approach can be used in a more general case of the “concep- 
tually stratihed” flow, see e.g. [BejanOS]. The main difference is that in 
this case an effective height for the liquid level is determined by equat- 
ing two expressions for the Lockhart-Martinelli [Lockhart49] parameter. 
The first expression relates Martinelli parameter to the geometrical char- 
acteristics of the flow (see [BejanOS]) 


Xtt 


^gwd T h 


IT 


^gwd T 
■^gwd 




A 


Iwd , 



1 

4 


43 

^Iwd 


42 

^gwd 


hwd 


(12) 


The second expression correlates the value of the Martinelli parameter 
with the thermodynamic characteristics of the flow as follows 


Xtt — 



(13) 


To solve the equations (12) and (13) with respect to bi the values for 
the phasic cross-sections are taken in the form [BejanOS] For bi < 0.5: 


Alwd — 


(8{biy^ + i2[bi{i-bi)f) 

15 


7T 


and Ag^(^ — 


For b[ > 0.5: 


(8{l - bi)-^ + 12[bi {I -bi)f) ^ 

Aguud and ^ Ag^^^. 

Here, all the geometrical parameters are made dimensionless (as in- 
dicated by additional subindex d) by scaling with the pipe diameter D 
or diameter square when appropriate. 

Once geometrical parameters are determined one can calculate source 
terms as discussed in more details in the following subsection. 


4.3 Heat fluxes 

As a first step in the analysis of the source terms we consider heat fluxes 
per unit volume at the dry and wetted wall and at the interface given 
by the following equations written for each control volume 
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Figure 9. Boiling curve example (water-vapor) for a given pressure, mass 
flux, and subcooling. 


An additional condition has to be imposed on these fluxes, which 
states that in the sum of energy equations the interface terms must sum 
to zero: 


+ ^il,L + ^7g,L{H?g,L ~ H^,l) ~ 0 - 


The full set of correlations required to calculate heat fluxes (14) for 
various flow patterns as a function of void fraction, wall superheat, and 
mass flux will be discussed in details in a separate report [LuchDG-III] . 
Here we provide a simplified discussion of the heat transfer correlations 
near the wall, neglecting heat fluxes at the interface. To follow the 
transition between various heat transfer correlations as a function of the 
wall superheat AT^g {AT^s = T^ — Tgat) let us recall the typical shape 
of the boiling curve shown in the Fig. 9. 

There are three characteristic temperatures that separates four re- 
gions with different physics of the heat transfer for AT^s > 0 (see 
e.g. [Bejan03,Nellis09]): 


Tonb > Tgat Convective heat transfer, which is characterized by com- 
plete contact of the fluid with the wall, and can be natural or forced, 
laminar or turbulent, single or two phase depending on the mass 
flow rate and mass fraction value; 

Tchf > Tuj > Tonb Nucleate boiling that occurs when the wall tempera- 
ture is above the temperature of onset of nucleate boiling (Tonb) 
and is characterized by bubbles nucleation, growth, and departure 
from the heated surface; 

Tmin > Tni > Tchf Transition boiling, which is an intermediate regime 
between the nucleate boiling and film boiling regimes that occurs 
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when the wall temperature is above the critical heat flux temper- 
ature iTchf)- The heat flux tends to decrease, while the dry wall 
area tends to increase with an increase of the superheat. 

^ Train Film bailing, in which a stable layer of vapor that forms be- 
tween the heated surface and the liquid, such that the bubbles form 
at the free interface and not at the wall. It occurs when the wall 
temperature is above of the Leidenfrost temperature Tmin- The 
heat flux tends to grow with the increase of superheat. 

The heat transfer coefficient is determined as 

h=^Nu, (15) 

where k is thermal conductivity of the fluid, D is the pipe diameter, and 
Nu is the Nusselt number. The correlations are given in terms of the 
Nusselt number. 

Convective heat transfer. In the single phase regions the follow- 
ing correlations for laminar and turbulent forced convection and natural 
convection [RELAP5-IV, TRACE] 

r 4.36, 

_ I 0.023 • 

“I 0.1 • {Gr • Prf/^, 

[ 0.59 -{Gr-Prf/^, 

are chosen for the Nusselt number for both gas and liquid. To guarantee 
a smooth transition between the various regimes the maximum of the 
above numbers is taken as the value for the convective heat transfer. 

Here, Re = ig Reynolds number for a given flow, Pr = 

and Gr = ^ 9 h(T,u Grashof number of the flow. 

The convective heat flux is finally determined as follows 

Qc = he {Tre - Ti) . 


forced (Lm) [Sellars56] ; 
forced (Tb) [Dittus30j; 
natural (Lm) [HolmanSOj; 
natural (Tb) [HolmanSOj. 


Nucleate boiling. To determine the heat flux corresponding to 
the nucleate boiling, one can follow e.g. [TRACE] and find first the 
temperature of the onset of the nucleate boiling Tonb 


Tonb ~ T ^ 


^ ^Tonb,s + ^Tonb,s + 


(17) 


where the correction factor E(</>) as a function of the contact angle 4>, 
subcool temperature ATgub, and saturated temperature of onset of nu- 
cleate boiling A.Tonh,sat are given by the following formulas 


AT sub = Tsat - Tl; ATonb,s = 


2hfcCrT 


sat 


F‘^{4>)PgHlgKf 


F{(f) = 1 - 
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Once the temperature corresponding to the onset of nucleate boiling 
under current flow conditions is determined the corresponding heat flux 
is found using the following simple correlation 


Qnb — 


Qc “1“ i^Qpb Qbi) 


1/3 


(18) 


where Qc = he iT^ — Ti) is the convective heat flux, found above, qp^ is 
the pool boiling heat transfer, qf,i = qpb{Tonb) is the pool boiling heat 
transfer at the onset on nucleation. To complete the calculations of the 
qnb we have to determine heat transfer coefficient for pool boiling using 
e.g. the following equation (see [TRACE]) 

1 71 

hpb = (ho • Fp/qo) i-" • (T^ - Tsat) , (19) 

where ho = 5600 [VE/m^/AT], qo = 2000 [W/m?], n = 0.9 — 0.3 • 

Fp = 1.73 • ^6.1 + ■ Pr^, Pr = P/Per, and Per is the 

pressure at critical point. 

Transition boiling. Transition boiling corresponds to the interme- 
diate regime between nucleate and him boiling. The transition boiling 
heat hux is usually given as a result of interpolation between character- 
istic heat huxes qehf and q-min- We will use the form of interpolation 
introduced in [TRACE] 


Qtb ftb ■ Qchf T (1 ftb'iQmin: 


( 20 ) 


where 



It can be seen from equation (20) that to hnd qtb one has to determine 
values of the four parameters: (i) qehf, (h) q-min, (hi) Tehf, and (iv) Tmin- 
There are a few different approaches to characterize these character- 
istic values. Eor example one could hnd the hrst two parameters using 
tables or correlations. Multiple correlations are available for the values 
of qehf and qmin [KandlikarOl]. One of the best known correlations for 
the critical heat hux was proposed by Kutateladze [Kutateladze48] 

Qehf = K • higpg . (21) 

A well known correlation for the minimum heat hux was suggested by 
Zuber [Zuber58] 

qmin = c • higpl/^ j . (22) 

Coefficients K and C in the equations (21) and (22) are of the order of 
0.1 for the water, but are not well known for the nitrogen [KandlikarOl, 
Yuan06] and can be used as htting parameters. 
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Once Qchf and qmin are found the characteristics temperatures Tchf 
and Train are calculated as follows 


Qnb{Tchf) — Qchfi Qfb{Tmin) — qmin- (^ 3 ) 


Alternatively, one could determine values of Tmin and Tchf using 
correlations and then apply equations (23) to find qchf and qmin- 

In practice, the set of correlations is optimized by comparison of the 
model predictions with experimental data as will be discussed in more 
details in a separate technical report [LuchDG-III]. 

Film Boiling. To hnd heat flux at the wetted perimeter the film 
boiling heat transfer is used in the form [BromleySO] 


hfb = C 


apgi^l {pi - Pg) HigCpg 

D {Tn, - Tspt) Pvg 


0.25 


where C = 0.62, TfUm = | {Tw + Tspt) is him temperature and 


Hlg =Hg- Hi 


is the effective heat of vaporization. It is assumed here that all the heat 
transferred from the wall to the liquid through the wetted perimeter is 
be used to heat it up and to evaporate. 

We note that the him boiling heat transfer is one of the key properties 
of the cryogenic how that affect the chilldown process. 

Flow corrections. 

For the how boiling the values of the characteristic heat hux q and 
temperature T have to be corrected taking into account mass hux and 
void fraction of the how as will be discussed in details in [LuchDG-III]. 


4.4 Mass fluxes 


The total mass transfer per unit volume in each L— th control volume at 
the n— th time step is dehned as a sum of interfacial mass transfer 
T” ^ and near wall mass transfer T”^ ^ 

■pn ^Y^n I -pn 

g,L — ig,L ~r wg,L- 


The interfacial mass transfer is given by the following relation 

€i,L + 


■pn 

~ - H 


g*,L 


l*,L 


where the denominator is given by 


Trn in _ 

^g*,L ~ — 


H^s,L-Hr,Li r>0; 
hIl-h?s,l^ r<0. 


The saturation values are determined using NIST tables [NIST90]. 
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The near wall mass transfer is given e.g. by [RELAP5-IV] 



max 



where correction coefficient M^or is discussed in more details in a separate 
chapter related to correlations. 

In the simplified version of the code we neglected the interface heat 


simplihcation can be only partially justihed by reasonable agreement 
with the experimental data during earlier chilldown stage and will be 
removed in the next version. As a further simplification accepted at the 
development stage the ambient heat transfer coefficients was taken as a 
constants over the whole temperature range. 

Once heat and mass fluxes are found one has to determine boundary 
conditions at the input/output valves of the system (including dump 
valves) as explained in wore details in the following section. 

4.5 Boundary Conditions 

The solution of the expanded equations (30) requires knowledge of the 
upwind values of the flow variables and velocities at the inlet and outlet 
of the system. 

Thermodynamic characteristics for liquid and gas in the storage and 
vehicle tanks are found using separate subroutines for each tank. 

Storage tank. Pressure in the storage tank is one of the main control 
parameters in the system and is considered as a given boundary condi- 
tion. For the storage the vapor phase is assumed to be at saturation 
temperature corresponding to a given pressure. During loading oper- 
ation the liquid in the storage tank is generally subcooled with liquid 
temperature being close to the equilibrium temperature at atmospheric 
pressure. 

Vehicle tank. The vehicle tank at the KSC testbed is ventilated at all 
time during loading operation. And there is no back flow to the transfer 
line from the vehicle tank. Accordingly, the boundary condition at the 
exit of the transfer line is determined by the atmospheric pressure and 
hydrostatic pressure of the liquid in the tank. 

To hnd gas and liquid velocities through the input and output valves 
one should use, in general, a two-phase flow model of the valve. Cur- 
rently, for the sake of simplicity the flow of each phase through the valve 
is assumed to be independent and incompressible, which is reasonable 
approximation for gas velocities less than 50 m/sec. The void fraction 
of this flow through the valve is assumed to be the same as the void 
fraction of the incoming fluid. The resulting volumetric flow rate is 


exchange in the bulk and considered only mass transfer at the wall. This 



(24) 
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Figure 10. Schematics of the pipe with dump valve. K 2 , and are 
flow coefficients for the dump valve, check valve, and other minor losses 
respectively. 


The coupling of the pipe flow to large volumes in the storage and vehicle 
tanks is modeled by taking into account the inertia of the flow through 
the input and output valves in the form 

Jvi = (25) 

'^vl 


where 


4.6 Dump valves model 

The mass flow through the dump valves is modeled using the following 
simplifying assumptions: the pressure at the inlet of the dump valve can 
be taken as a pressure in the control volume coupled to this vale, while 
the pressure at the outlet of the dump valve is approximated by the 
atmospheric pressure. These approximations are justified by the short 
pipes of low resistance connecting dump vales to the transfer line and to 
the drain system. 

The mass of the gas flow through the valves can be approximated in 
two different ways. In one of the approximations the flow is considered 
to be compressible. As a result the following equations for the mass flow 
rate can be used for the dump valve 


= S.i 





r, .2 

7 + 1 

27Pin 

/ Pout \ T 

/ Pout \ 7 

7-1 

V Pin ) 

V Pin J 


in supersonic regime 
in subsonic regime 


For relatively low gas velocities the incompressible approximation of 
the flow through the valve considered in the was found to be accurate 
enough for practically all the loading conditions at the KSC testbed. 
Typical conhguration of the pipe with dump valve is shown in the Fig. 10. 
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The corresponding volumetric flow rate through the dump valve can be 
found as follows 


Qo 


f I 11 

+ iq + Ki 


- 1/2 



PLN2 

P 


where ci is the relative opening of the dump valve. 

To take into account the inertia of the valve operation characterized 
by the time delay Ty the volumetric flow rate through the dump valve 
was modeled in the following form 


Qo — Q 

TV 


The heat flux through the dump valve H^v was then calculated as 


Hdv — Q Pg^p^gi 

where CpTg is the gas enthalpy in the control volume attached to the 
dump valve. 

Once the boundary conditions including heat and mass fluxes at the 
liquid/ gas/ wall interfaces and through the input /output valves are found 
the algorithm proceeds to the calculations of the first integration step. 


5 First step 

The first integration step of the algorithm includes the following sub- 
steps: 

1. Solve expanded equation with respect to pressure in terms of new 
velocities; 

2. Substitute this solution into momenta equations and solve resulting 
block tri-diagonal matrix equation for the new velocities; 

3. Find new pressure; 

4. Find provisional values for the energies and void fractions using 
expanded equations; 

5. Find provisional values for mass and heat fluxes. 

We now discuss in more details each of the sub-steps. 

5.1 Expanded equations 

The set of expanded conservation equations discretized on the main grid 
(see Sec. 2.2) include: 
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the sum density equation 


+ KMX" + ^< 1 ' {pIl - pIl) + 
{(WA)lt^M»-KpA)l,uM + 

{fppA)M Kxx - kpA)i, = 0 , 


At 

V 


At 

V 


(26) 


the difference density equation 


p;.MX" - KMX" + (pSa + p"l) + 
f ((apA)rj^,urXX - (SpA)rjKX") = 


At 

V 


2r 


ff,L’ 


(27) 


the gas energy equation 


(Pg,L®g,L + Pl) + (^P)g,L '^Pg,L^ + 


y 


(Aa (pe + p))gj+i M^++i + (^a (pe + p))”^- u' 


72+1 


H. 


{ rpn ^n+l\ Q I TTn ( rJ^s,n+l ^n+l\ q 

gw,L \^w,L ^g,L ) ^wg K- J^ig^L ^ g,L ) 


_l_-nn jAn _i_ -pn jAn 
ig,L^ig,L “I" wg,L^wg,L 


At 
V ’ 


(28) 


the liquid energy equation 


- (+L+L + Pl) da^^L + MIl + («e)+ dpi 


,n+l I 
L 


At 

V 


(Aa {pe + p))lj+i + {Aa {pe + p))lj u\ 


n+l 


jAn 
^lw,L 


{tS.l - KX") Ki + ifS,i (TiX+' - r,y‘) s, 


pn trn ■pn jAn 
^ ig,L^il,L ^ wg,L^wl,L 


il,L 
At 
V 


(29) 


5 . 1.1 Solving expanded equations with respect to pressure 

The four coupled equations (26) - (29) for the sum and difference density 
and for the gas and liquid energy can be rewritten in matrix form 




dcxg 

deg 

dei 

dp 




= a 


MX"+y + l>MX"+‘^dulp, + didufp + el, ( 30 ) 


Here the vector of unknowns = |iiQ;”+, cie”+, 

Assuming for now old values of the temperature on the right hand 
side of eqs. (26) - (29) and using the following expansion for the densities 
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the elements of the matrix can be written in the following form (index 
(j) enumerates columns of the matrix A^) 


= 


1 

ej • 

1 


~^Il {deP)lL 

«g,L {iodep)'(,,L + P'^,l) 

II 

SIh 

0 

0 

/3g,L [ppL + iedep)l,L) 

«g,L (^eP)g,L 


^II ideP)lL 


pIl + P?,l 


■ <L {dpP)lL - {9 pP)Il ' 

{opTg^L+Pl 

II 

{eadppfgj^ 

- {op)Il - Pl 

{fdedpp)'^ ^ 

. PIl - P?,L . 


. {9vP)Il + {9pP)Il _ 


= 

■Ax ry^ 


Columns of the vector-multipliers for gas velocities on the right hand 
side of the eq. (30) have the form 


< = - 




(ep) +pI 


ff J+1 
0 

Pg,j+^ 
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'9J+1’ 
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(ep) +pI 

V / 0.7 


9,J 
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.^n 

Pg,3 
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g,j' 


Vector-columns for the liquid velocities are 



r -—n n 

Pl,j+l 


r -I 

~ Pid 


0 


0 

= 

^x 

- (ep) ~Pl 

. - r!^ — 

(e'pY +pI 


-Pid+i 


L Pid -1 




where the coefficients nf,- are of the form 


^tA 

^9d = 


3 . 


nij = oil, 


AtAi 


-P3 Vl ■ 

Finally, the free vector in eq. (30) is written as follows 


e” = 


^llHZgVL - Q^g 
-^IlHZiVl - Qn.1 
0 


At 


' ^x^gj ' ^x^lJ-\-l ' ^x^lj 


The matrix equation (30) is solved with respect to dp^^ in terms of 
new velocities dUgj^i, du^j^, du^~^^ . 

The new pressure dp^^ is expressed in terms of new velocities 

du!l~j^ using the last row of eq. (30) multiplied by the 
inverted matrix A” 


dp 


n+l 

L ~ 




V7/”+^ 
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where is the 4-th row of the inverted matrix Ax- 

We note that excursive numerical tests has revealed the high sensi- 
tivity of the cote to the accuracy of numerical inversion of the matrix 
Ax. This is due to stiffness of the problem. Indeed, the eigavalues of the 
matrix are different by twelve orders of the magnitude. As a result 
the errors of numerical inversion of matrix A” cause significant stability 
and accuracy problem. To solve this problem we used the exact equa- 
tions for the elements of the matrix in the code. Such a simple 

solution sufficed to improve the stability and accuracy and also allowed 
to speed up the calculations. 

To build matrix A” and vectors a”, 6”, c”, and e” ones has to 
calculate 

• upwind values for void fraction, energy, and density, 

• density derivatives with respect to pressure and energy, 

• source terms, i.e. right hand sides of eqs. (26)-(29). 

The analysis of the heat and mass fluxes required for the calculations 
of the source terms was given in the Sections 4.3 and 4.4. Below we define 
the upwind values and discuss calculations of the density derivatives. 
Note, that in the current version of the code the temperatures in the 
right hand sides of the energy equations (28) and (29) are taken at the 
previous time step to simplify the calculations. 


5.1.2 Upwind values 

The values of the velocities, mass, energy, and void fractions are cal- 
culated at the interfaces centered at staggered grid with indexes j = 
1, . . . , -|- 1. The values of scalar variables are calculated using upwind 

scheme following [RELAP5-I]. For the density we use 


= 5 i^L-l +y^l) + l {Su,j + ZuJ ■ Spj) {(fL-l - ^l) = 

(1 T (^uj T ZuJ ■ ^pj)) 2 ^ ) 

For the energy and void fraction we use 

= 2 (^L-l + V'l) + 5 {Su,j + Zu,j ■ Spj) {ipL-l - V'l) + 

Zp^j 2 {tpL-l + VL) + P,(^i^,L-l+Pg(l),L J ~ 

(l -h {Su,j + Zu,j ■ Spj) - ZuJ ■ ZpJ + 

^1 T ZuJ ■ Spj 


Sr).j ) Zn, 


‘“J + PgJ),L-l+Pg(l) 

Here the following notations are introduced 


2Pg(0,Z,-l 


A 


2 • 


(32) 


( 33 ) 


Suj = sign{uj)] Zuj = {1 if uj = 0, 0 otherwise} ; 
Spj = sign{pj); Zpj = {1 if pj = 0, 0 otherwise} . 
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T, K T, K 


Figure 11. Comparison of the values for the density derivatives obtained 
using ideal gas equations (red circles) and [BridgmanOl] equations (blue 
dots). 


It is also recommended (see [RELAP5-I]) to use additional smoothing 
for small mixture velocity jm = ^gjUgj + aijuij in the form 

{ 1 when jm > Jiim 

(3 - 2^) when - jum < jm < Aim 
0 when jm < -jum 

where ^ and Aim = 0.46 m/s. In the present version of the 

code this smoothing was omitted to simplify the calculations. 


5.1.3 Density derivatives 

In the present version of the code the temperatures in the source terms 
are taken at old time. Therefore, only density derivatives have to be 
calculated. Density derivatives are found using two-dimensional interpo- 
lation of NIST tables [NIST90]. The required values of the g^[T] could 
be taken directly from the tables. In more general case the derivatives 
were estimated using the following equations [BridgmanGl] 

dejp (^cp- \dp)u (cp - Pp/p) 

In addition, at the development stage the density derivatives for the 
gas phase were frequently approximated using ideal gas equations 

(f) =r^- 

\dejp \eg J \opJjj (7-l)eg 

The use of ideal gas approximation can be justihed by direct com- 
parison of the derivatives values obtained using [BridgmanGl] equations 
and ideal gas equations as shown in the Fig. 11. It can be seen that ideal 
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Figure 12. Comparison of the values for the density derivatives obtained 
using table data (red circles) and polynomial approximation (black solid 
lines). 


gas equations approximate table values with reasonable accuracy all the 
way down to the temperatures around 100 K. 

For density derivatives in the liquid phase simple polynomial approx- 
imations can be used as shown in the Fig. 12 

We note that for any type of approximation a special care has to be 
taken to avoid singularities at critical temperature. 


5.2 Momentum equations 


One of the critical steps of the nearly-implicit algorithm that determines 
its accuracy and speed is the solution of the velocity matrix, which we 
will discuss now in more details. 

The original non-conservative version of the gas momentum conser- 
vation equation for the Wallis model has the following form (see (1) and 
(2), cf [Staedtke], [Wallis69]) for the gas 


AapgUg^x + \Aapg[uf^ ^ + Aap^x = -AapgZ^x- 

(^) - CiUR |njj| + ATg (ugi - ug) + Mv 

and for the liquid momentum conservation equation 

^l3piUl,x + l^PPl{u‘l)^x + = -APPl^,x- 

- aUR \ ur\ - ATg {Uli - Ul) - My. 
Here My is the virtual mass term 


My 


Ca (1 - a) Pm 


d{Ug - Ul) 

dt 




dui 
^ dx 


(34) 


(35) 


( 36 ) 


and Pm is the mixture density 


Pm — ^gPg ^lPl‘ 
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5.2.1 Numerically convenient form 

To enhance the stability of integration in semi-implicit and nearly-implicit 
algorithms we follow recommendations of the RELAP5 code and use in- 
stead of the original equations (34) and (35) the sum 


aPgUg^t+ iSpiUl^t + ^ = -PmZ,x- 

aPgUgFyjg ~ / 3 piUlF^l ~ Tg {Ug ~ Ul) , 

and difference momenta equations 

(v?„) (uf) / \ 

Un,t ~ Ul^t H \Fg Fi) ~ ~ UgFyjgF 

^ {Ui - Ul) + PmFi {Ug - Ul) + ap^Ppi 


apg ^ 9 ) + Bo, ('“'9 + aoTeoi^'^ 


(37) 


(38) 


The following simplifying notations were introduced for the wall fric- 
tion (cf [RELAP5-I]) 


Aoig PgUgFugg fglwg 


PgUg 


Un 


AaipiuiF^i = filwi 


PlUl \Ul\ 


(39) 


Similarly for the interfacial friction we have 


AapgFig {ug ~ Ul) = A {l ~ o) piUiFii {ug - Ul) = CgiUR \ur 


Note that in writing sum and difference momenta equations the orig- 
inal equations were divided by cross-sectional area A. In addition, in 
deriving the difference momenta equation the original equations were di- 
vided by the product of void fraction and density for each phase. Note 
also, that difference equation involves reduced virtual mass term My, 
which excludes spatial derivatives in ( 36). 

In the nearly-implicit algorithm the kinetic energy terms are treated 
approximately implicitly as follows 



n+1 

g,L 


u 


n 

g,L 




+ 



2 


The resulting finite-difference momentum equations for one control 
volume on the staggered grid can now be re-written for the sum (40) 


' \"rn,3 (j 

f {ap)l. + u^gl - 2au-F-iUg,L ~ % 

- {dpt^ - dplt\) At - [pl - Pl_,) At- 


- ‘^du^gt-lUlL - <;Ll) 


AtAxj 


PmjPAz^ + Mlj Fn,g,g ^ ^^g,j ) ^ 

.i K/‘ + “u) - Kj‘ + "h - ‘'“"b - <i) 


n jpn 


.9 


+ <A + 


wl 


( 40 ) 
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and for the difference (41) 


(i + si)) Ki' - ‘‘“ir] 

^ ((i)Ij + “is “ - “is-i 

- ( W) (fi*‘ - 

/ \ P" .d<+l-(op)'* .d<+l-(ap)".d«'*+ 

Tpn I ^,.n+l I ^.n ^ tti ^m,j i,j i ^^g,j l,j '■ ^^1,3 9,3 

^Wlj [d^l j + Ul j) - i gj (op)", (op)". 


(ap)",-(«P);;,- 

(F,rt” (<!<+' + - <i<+‘ - ufj) } AtAx, 


n+1 


(41) 


,3_gjj_ ^ 


of momenta equations. 

The matrix form of the last two equations reads 


b: 


du 

dv 


“I n+1 


+ CJi 


Jj-i 


du 

dv 


“I n+1 


+ D'' 

I -L^ 7 , 


du 

dv 


bo + bidpl^^ + b2dplt[ 


-‘j 


n+1 


-I n+1 


Ji+1 


(42) 


It can be noticed that the relations between the new velocities is 
non-local for the two main reasons. Firstly, it involves velocities at faces 
{j — l,j,j + 1}. Secondly, it depends on new volume centered velocities 
Ug^k and momenta pk at locations {L, L — 1}. 

In turn, volume centered pressures and velocities are related to the 
face centered velocities at the neighboring cells. This relation for the 
pressures was obtained earlier (see eq. (31)) and is repeated here in a 
more compact form 

dpl^^ = o^pdUg^h + bpdu^+^ + c^du^+l^ + d^du^^^ + e^, (43) 

where coefficients on the right hand side are Op = a”, bp = 

(T-)J' 6-, Cp = (x4-)ji +, dp = (x4-)-i d-, and Cp = (A-)^' 

The corresponding relation for the volume centered velocities du^'^^j^ 
in terms of the velocities at the faces is given in the next section. 


5.2.2 The volume-centered velocities 

The volume-centered velocities can be in general represented in the form 


'^g{l),L ^g(l),L'^g(l),j+l + bg(l),L'^g{l),j (^4) 

For example for a straight section of the pipe we have 


%(0T 


I ++l®g(Z),i+lPg(Z),j+l _ 

o A ^gd)J ' o A ^p(0 J+i 

^^Lag^l),LPg(l),L 2ALag^i)^LPg(l),L 


(45) 
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5.2.3 The velocity matrix 


On substituting (43) and (45) into (42) and combining the latter equa- 
tions for all N faces we obtains the velocity matrix equation in the fol- 
lowing form 
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C 12 
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^21 
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„N-1 
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Oil 
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dug^ 
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d^^N 



Ui 


ri2 


ng 
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^2N-2 

"2^-1 

^2AT 


5.2.4 The block tridiagonal solver 

It can be seen that the velocity matrix has the block-tridiagonal structure 
of the form 


Bi Cl 
A 2 B 2 C 2 

■^N-l Bm-1 Cn-1 
An Bn 



Xi 


ki 


X2 


k2 


XN-1 


kN-i 


XN 


kN 


(46) 


Here A, B, and C are 2x2 matrices, Xm = 

iT 




-\T 


and km = 


n. 


n 


2m- 1) ^2m 


n 


are 2-dimensional vectors. 


Applying equations (43) and (44) to momenta equations (40) and 
(41) we obtain the explicit form of the matrix coefficients in (46) given 
below. For the coefficients of the Ai matrix 


«Mi = bg,L-iulL-i - Op2,L-i) At; 

«fc,12 = h,L-lUlL-l - Op4,L-l) At 


(47) 
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«fc ,22 = 11 ^ 




For the coefficients of the B{ matrix we have 


un _ 


_ (' + + r;,i) Al;+ 

+ (“p2,L - Opl,L-l)) At; 


°fc ,12 — 


(((^)i:,(i + i^,:j/,,)-rL) Ax. 


( 48 ) 


(ap)”^. (hig^LulL - ai,L-iulL-^ + {ap^,L ~ flpS.L-i)) At; 


/.n _ 

"fc ,21 — 


(l + )jAxj + At {bg^LUg^L ai^L-iuf^L-i) 


+F" , Ax, At - 7=^ - 1 Ax, At 


wg,j 
+ 


(ap) 


A"p" 

T " rj 




(^I^) («p2,L - api,L_i) At + p^jF^jAxjAt-, 

(l + ‘®^),Axj - At [KluIl - ai,L-iulL-i) 


— — 
^k ,22 — 


TT'^ /\ T* • /\f ^ 

{ap)l. 


^±^^-'-SAxAt+ 


ap 3 ,L-i)At 

For the coefficients of Cj matrices we have 


. /\ T • /\ f ■ 


4,11 = U(^PTg,j «9,i^g,L + Ophi) At; 

Cfc,12 = {i^P)lj + «p3,l) At; 


''fc ,21 — I I — 


^k,22 — I I — 


{^)gj^9FUg,L + ^phLj At; 

(ir,«^^<L+(»)«p 3 ,L) At 


(49) 


The free vector has the form 


"2t = *j9 (%f ) Aj/2Af - ((«;.i)" - (t.;,i_i)") 

^((“”l)'-(“"i-i)') 
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{pl - Pl_, + At - (u; 
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Many methods exist that solve block tridiagonal equations of the 
form (46). We are specifically interested in the extension of the Thomas 
algorithm [Reimar66], which allows to use the fact that our matrices 
have size (2 x 2) and inversion and multiplication of the matrices can be 
explicitly coded. In [Reimar66] it is demonstrated that a simple forward 
elimination-backward substitution is successful (neglecting rounding er- 
rors for the moment) if all minors of the matrix in (46) are non-singular. 

The algorithm consists of two steps. Forward elimination 


Hi = -B^^Ci 

Hn = —[Bn + AnHn-l] Cn, n = 2, . . . , N — 1 

gi = B^^ki 

9n — [Bn T -^nHn—ll {kn -'^nOn—l') ; 7T, = 2, . . . , N 


(51) 


Backward substitution 


XN = 9n 

Xn — 9n T HnXn+l) U, — N 


(52) 


In these notations the explicit form of the matrix coefficients 


A 


n 
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- n 

a\ as 

02 U4 

J j-i 


Bl 


bi 

b2 



Cl 

C2 


cs 

C4 


1 ^ 


Jj-hl 


where matrix coefficients are given by equations (47) - (50). 

To find the solution for the new velocities the following variables in 
the equations (47) - (50) have to be calculated 


• face centered values of the thermodynamic variables; 


• coefficients of the virtual mass; 

• friction terms; 

• velocities at the input/output valves. 

The corresponding calculations are discussed in the following sec- 
tions. 


5.2.5 The face-centered values 

The face-centered thermodynamic variables are interpolated between 
neighboring cell values based on the length of each half cell 

+ 

where represents densities and void fractions for liquid and gas. 
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5.2.6 The coefficient of the virtual mass 


Following [RELAP5-I] we introduce the coefficient of the virtual mass in 
the form 


Mv = 


1 (l+2a) 

2 1-0 

1 3-2o 

2 a 


for 0 < a < 1/2, 
for 1/2 < a < 1. 


The common multiplier for this coefficient remains a fitting parame- 
ter established on the basis of numerical experiments. 


5.2.7 Friction terms 

One of the difficulties of the pressure based solvers on the staggered grid 
is the approximation of the momentum equation for abruptly changing 
pipe diameter. A quasi-steady approximation is used in [RELAP5-I] 
and continuity of the volumetric flow is used in [TRACE] to mitigate 
this problem. 

In the current simplified version of the code the following approach 
is adopted. The friction coefficients fg, fi, and fi in equations (37) and 
(38) are determined using Churchill [Churchill77] approximation for all 
flow regimes (54) 


U = 2 



1 

(a + 6)3/2 


-1 1/12 


(54) 


with Reynolds numbers 




Pm,LUm,LD 

m^L 


based on volume centered velocities Um,L given by (45) and hydraulic 
diameter 


D 


m 


4 Al 

^m,L 


for each control volume. Here m takes values m = {g, I, i} for gas, 
liquid, and interface in a given control volume. 

Next, the major losses are presented in the form (39). And the minor 
losses are given by 





2 


Here we have assumed a simple version of partitioning minor losses be- 
tween phases. According to this equation for the sum of the momenta 
with equal velocities the homogeneous mixture has the same losses as 
liquid single phase. A more elaborated versions of the minor losses par- 
titioning in the two-phase flow will be considered later. 
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With the simplified partitioning the sum of the major and minor 
losses for a given control volume can be written in the form 

2 I I 

{fmlm^X + aniAKml) + K^i) A OmPmUm-^ , (55) 

where 

r-i n |^m| 1 j, f-ml 

Xm — Jm T and -Uml — Jm~pr ■ 

amA 1 D 

Here Imi effective pipe length corresponding to minor losses. 

In the proposed above simplihed approach minor losses are added to 
the major losses. Partitioning of the minor losses in a more general case 
will be discussed elsewhere. 


5.2.8 Boundary conditions for the velocity matrix 


Boundary conditions for the velocity matrix are determined by the liq- 
uid/gas flow through the input/output valves. The flow through the 
valves is assumed incompressible as was discussed in Sec. 4.5 and the 
corresponding 2x2 boundary matrices A, B, and C in the velocity 
matrix are of the form 


A = 


0 0 \ 

0 0 j ’ 


/ \ 

1 ^ ^ ’ 


C 


0 0 \ 

0 0 j ’ 


while the vector k is written as follows 


k 


„,n+l I „ n+1 K+npdt 

“T ^10 T 

,,n+l n+1 (u^-u")dt 

“gO ^10 T 


where 


n+1 

^9(00 



/^pU+l 



(56) 


where A^ is the pipe diameter fitted to the valve, Ky is the flow coeffi- 
cient, Ap is the pressure drop across the valve, and Pg{l) is the density of 
the fluid flowing through the valve. It is assumed that the void fraction 
and the temperature of the fluid are the same as in the upwind control 
volume. 


5.2.9 Structure of the solver for velocities 

To complete this section we briefly outline the structure of the solver for 
the new velocities. The following values are calculated at this sub-step: 

1. Face centered values of thermodynamic variables; 

2. Source terms for friction and inter-phase mass fluxes 

• Gas wall friction; 
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• Liquid wall friction; 

• Interface friction and 

• Face centered interphase mass fluxes. 

3. Coefficient of the virtual mass; 

4. Face centered mixture density; 

5. Face centered density difference; 

6. Volume centered velocities; 

7. Boundary conditions; 

8. Coefficients for matrices in block diagonal matrix Ai, Bi, and C*; 

9. Coefficients of the free vector; 

Once these calculations are completed one can solve block tridiagonal 
matrix and update velocities. When new velocities and new pressures 
are found, the provisional values of the densities, energies, void fractious, 
and temperatures can be found as described below. 


5.3 Provisional values of thermodynamic variables 

To find provisional values of thermodynamic variables new velocities are 
substituted back into equation (30), which is then solved with respect 
to Solution for the pressure is 


obtained in the form 

= Op + bpdu^j^ + c^du 

and for the void fraction in the form 


n+1 

tg+i 


+ dpdu 


n+1 

bj 


+ e 


L 
P ’ 


With a, =_(T^)r'a", ba = c, = c^, da = 


e, = (A-)-i ei 

Similarly, we have for the change of the gas energy de^'^ 

= ^^gd'^gj+i + + ^e,gdu"ljli + d^^gdu!l^^ + e^^g, ( 58 ) 

with ae,g = (41-)-' a-, 6e,g = (A-)^-' 6- Ce,g = (41-)"' C- de,g = d^, 

ee,3 = (4LS)2-'e- 

and for the change of the liquid energy de^~j^^ 

de^X^ = ae,idUg^h + b^,idv^f + (X^idu'ljli + d^^idu^X^ + ( 59 ) 

with ae,z = (41-)-' a-, = (4l-)3-' 6-, Ce,z = (41-)-' C-, de,i = (41-)-' d-, 

ee,/ = (A-)-'e- 
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Once values of the energy and pressure are found one can calculate 
provisional values of densities and temperatures. For the densities we 
use the following expressions 



similarly for the temperatures we obtain 



Using provisional values of the densities and temperatures one can 
update the values of the heat and mass fluxes following the results of 
Sections 4.3 and 4.4. The updated intermediate values of the mass and 
heat fluxes can be used in the second step of the code. The details of 
the calculations at the 2nd step are provided below. 


6 Second step 


At the second step of the algorithm the unexpanded form of the gas 
density 
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and liquid energy 


d{ape)lt^ + (d (ape)”++i + (ape)"^+i + (ap)”j+i 
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(67) 
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equations are solved implicitly using values for the new velocities and 
pressures, and provisional values for the temperature, heat transfer co- 
efficients, and mass fluxes obtained at the previous step. 


6.1 Tri-diagonal form of the matrix equations 

The efficient solution of the equations (64) - (67) can be obtained by 
noticing that all four equations are independent and have the same struc- 
ture 

dUl+^ + = Sl, (68) 

where and the upwind values of the conservative 

variables can be represented in the form 

dLT+i = X-^+^dUl+l + n]+^dU2'^y (69) 

The eq. (69) for the density and for the energy has the form of equation 
(32) and equation (33) correspondingly. 

By substituting (69) into (68) all four equations (64) - (67) can be 
re-written in the standard tridiagonal form 

aLdU^ll + hLdU2+^ + CLdUlXl = Sl, (70) 

where = 1 + ~ = 

coefficients of the sub-diagonal, main diagonal, and super- 
diagonal of the tridiagonal matrix. 

The structure of the Sl terms is clear from the equations (64) - (67). 
Here we provide as an example the structure of the Sl term for the first 
equation (64) 


Sg,L = 


Vl 


- {otp) 


n n+1 

9J+1 '9J+1 


+ {ap) 


■ 


(71) 


The tridiagonal equation (70) can be efficiently solved using standard 
tridiagonal solvers (see e.g. [Thomas49, LuchDG-I]). Note that coeffi- 
cients aL^ bL, and cl for the tridiagonal matrix are the same for the 
density and energy equations. These coefficients are calculated only once 
for each pair of the conservation equations. 

Once conservative variables are calculated the density and energy for 
the gas and for the liquid can be found as follows. 


6.2 Calculation of the main dynamical variables 

Obtained conservative variables for the density and energy 

{MTl ^ (ape)”+^} 
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are used to find primitive variables 
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Note that in the fist equation above the provisional value of the liquid 
density is used to obtain liquid void fraction. 

It is important to emphasize that it is the efficiency of the three main 
solvers (for the set of expanded equations (26) - (29), for the velocity 
matrix (46), and for the set of unexpanded equations (64) - (67)) that 
render the nearly implicit algorithm as one of the most efficient methods 
for solution of the two-phase flow problem. 

For the one dimensional two-phase flow problems the number of it- 
erations scales linearly in time for this algorithm. For cryogenic loading 
system with up to 100 control volumes the integration time is a few 
hundreds times faster than the real time. Short integration time paves 
the way to the on-line optimization of the mitigation strategies making 
nearly implicit algorithm very attractive for applications to the on-line 
control of loading operations. 

However, stability of the algorithm poses a long standing problem 
of the two-phase flow research (see e.g. [Nourgaliev, LuchDG-I]) and re- 
quires a special attention. Stabilization techniques that are employed 
in the current version of the code are discussed briefly in the following 
section. 


7 Control 

The control module consists of three main blocks. In the first block the 
dynamical variables are reset to the limiting values when these values are 
outside of the predefined range. In the second block the values of the dy- 
namical variables are smoothed during phase appearance/disappearance. 
And in the third block the time step is reset when the deviation from the 
mass conservation is detected and when the values of the thermodynamic 
variables are changed abruptly or lie outside the physical range. 

In addition, the opening of the dump valves is controlled to limit the 
mass losses during one time step by values less or equal to 10% of the 
remaining mass. 

7.1 Direct control of the dynamical variables 

At the first control step the values of the thermodynamic variables are 
checked against their saturation values. If the phasic density or energy is 
beyond corresponding limiting saturation value it is reset to the satura- 
tion value. This simplification is adopted because in the present version 
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of the code the subcooled gas state and superheated liqud state are not 
allowed. In the future versions of the code the correlations module will 
be extended by including superheated liquid and sub-cooled gas states. 

At the first control step we also check void fractions against the 
minimum and maximum values. If the void fractions are beyond the 
limiting values (that are usually set at 1 x 10“^ and 1 — 1 x 10“^) they 
are reset to the corresponding limiting values. 

7.2 Smoothing of the dynamical variables during phase 
appearance/disappearance 

At the second step of the control module the void fractions are checked 
against the smoothing margins. In this work, we follow recommendations 
by Liou [Liou07] and adjust temperature, velocity, and density according 
to the following expression 

4>adj = g{x)4>d + (1 - g{x)) (72) 


where 

g{x) = (2x — 3) ; and x = — 

Xmsbx S-min 

Here “d” stands for disappearing phase and “c” for conducting phase. 
The role of the limiters and smoothers described in this section is crit- 
ical for the code stabilization. The exact values of the minimum and 
maximum void fraction, for which smoothing (72) is applied should be 
established using extensive numerical experimentation. Currently these 
values are set at the I x 10“'^ and 1 x 10“^. 

7.3 Time step control 

At the third step of the control module a time step of the algorithm is 
reset. Multiple checks are performed that can through a flag to reduce 
the time step of the integration. Firstly, it is insured that the changes of 
the densities, energies, and pressures during one integration step do not 
exceed 25% of their absolute values at the previous time step. Secondly, 
we test if the results of the calculations for the pressure and temperature 
are outside of the predefined range of values. Thirdly, we check if the 
predictions of the density and the temperature obtained using the Taylor 
expansion during the 1st step of integration do not deviate too much from 
the corresponding corrected values obtained by solving unexpanded set 
of conservation equations during the 2nd step. Finally, we control the 
deviation from the mass conservation by monitoring the overall changes 
of the mass in the transfer line and the total mass flow in and out of the 
transfer lime through the valves. 

For example, let us discuss in more details the control of the deviation 
of the corrected densities from the predicted values. The total density 
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for each control volume 



— a 


n+l 

9,L 



is calculated using solution of the unexpanded form of the mass and 
energy conservation equations at the second step of the algorithm as 
shown in Sec. 6. These total densities are compared to the densities 
obtained by Taylor expansion 



The maximum error for each control volume and the total error are found 
as follows 


= max 
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Pm,L 
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rm,L 
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Pm,L 
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Et(atd' 


If either of these errors is larger than 10 ^ the time step is halved and 
the integration step is repeated. 


7.4 Dump valve control 

Since the flow through the transfer line is actively controlled by the dump 
valves we have to ensure that there is no liquid flow through the damp 
valves in nominal loading regime. This is done in two ways. First, if the 
gas void fraction in a given control volume is between 0.2 and 0.05 the 
flow through the dump valve is smoothly reduced to 0 using smoother 
(72). Next, we check if the mass loss through the dump valve during one 
time step is smaller than 10% of the mass of the gas in the associated 
control volume. If it is not the time step is halved and the integration 
step is repeated. 


7.5 Code termination 

During integration the values of the key thermodynamic variables may 
occasionally be found to lie outside predefined range as described above. 
If this happens the time step is halved. The reduction of the time step 
continues until reached the minimum predefined value (usually 10 mks). 
If all the field values are found to be inside the range, and both mass 
errors are smaller than 10“^, the time step is reset to the maximum 
predefined value. The time step of integration is usually between 10 ms 
and 20 ms. 

If the integration time is smaller than a predefined limiting value 
and one of the errors discussed above still prevails the integration is 
terminated. 
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Extensive numerical testing shows that the control of the thermody- 
namic values and of the integration time step outlined above guarantees 
a stable and robust performance of the algorithm during cryogenic load- 
ing in a wide range of the system parameters with time step up to 20 
ms. Despite the simplified set of correlations the results of integration 
are in reasonable agreement with the experiment data. The comparison 
with the experimental data will be discussed in more details in the next 
section. 


8 Comparison with the experimental data 

In this section we briefly describe some of the key features of the al- 
gorithm and discuss a comparison of the model predictions with the 
experimental data of chilldown in the cryogenic transfer line. 

8.1 Some key features of the algorithm 

Partial verification and validation of the nearly implicit algorithm was 
discussed in our earlier report [LuchDG-II]. The original simplified ver- 
sion was written in MATLAB. Its main goal was to prove the capability 
of the method to deal with cryogenic loading and chilldown problems. 
The integration was performed on a uniform grid and the integration 
time was of the order of 80% of real time. 

Using this simplified version we demonstrated the capability to simu- 
late pressure waves and material waves, in particular, gravitational waves 
in the pipes. We demonstrated the phase separation due to gravity 
force, emptying the pipe filled initially with liquid, and vapor locking 
phenomenon in the regime when initially empty pipe is filled by the 
evaporating liquid. We also demonstrated the code capability to simu- 
late valves closing and opening and mitigation of the vapor lock fault 
during chilldown by opening dump valve. Finally, we performed initial 
validation of the model by comparing the results of simulations with 
the experimental data obtained during first 300 sec of chilldown of the 
cryogenic testbed. 

The C-|— |- version of the algorithm discussed in the present report 
preserves the properties of the original code listed above and offers a 
number of improvements. We have added the ability to integrate two- 
phase flow equations on an arbitrary non-uniform grid. The code was 
restructured and optimized to allow for much faster integration. The 
upgraded version currently runs up to 400 times faster than the real 
time. Such an acceleration paves the way to the on-line optimization 
of the mitigation strategies, which is important for efficient autonomous 
control. 

Another feature of the accelerated code is its ability to perform ef- 
ficiently search in the multi-dimensional parameter space of the corre- 
lation functions. This is essential for autonomous control of loading. 


40 


because the corresponding correlations are not well-known for cryogenic 
fluids and continuous on-line learning of the correlations is required to 
improve accuracy and reliability of the fault detection, identification and 
recovery. 

The analysis of the full set of correlations and their effect on the accu- 
racy of the model predictions will be given in a separate report [LuchDG- 
III]. A brief preview of required correlations that include flow patterns 
recognition and calculations of the frictional losses, heat and mass trans- 
fer was given in this report in Sections 4 and 5.2.7. Below we describe 
briefly a comparison of the predictions of the accelerated model with 
experimental results using a simplified set of correlations. Despite the 
adopted simplihcations the agreement with the experimental data was 
substantially improved as compared to the results reported [LuchDG-II] 
and extended from 300 sec to the 1600 sec of chilldown. 


8.2 Loading regime 

During chilldown at the KSC testbed the pipes and the vehicle tank are 
initially filled with hot nitrogen gas at T = SOOi^T and p = 1 atm. The 
storage tank is filled with liquid nitrogen at T = 80K and p = 3.245 
atm. At time 0 the input valve MV151 is opened manually. It takes 
approximately 15 sec to open this valve, however, the exact opening 
dynamics is unknown. The main control valve Roll5 remains closed for 
another 195 sec preventing flow through transfer line. During this pre- 
chilldown time two dump valves (dcvll2 and dcvll7, see Fig. 13) are 
opened at approximately at 145 and 163 sec to chill the first section of 
the transfer line. 

At approximately 195 sec the main control valve Roll5 is opened and 
the line is chilled and filled with liquid nitrogen during next 1400 sec. 
The flow variations depends mainly on the boiling dynamics and on the 
relative opening of the dump valves. The flow boiling dynamics is highly 
non-trivial process and will be considered in details in a separate report 
[LuchDG-IIIj. The valve opening, on the other hand is well defined and 
controlled remotely by the operators in the control room. Therefore an 
accurate simulation of the valve opening is a prerequisite for inferring the 
parameters of the flow boiling in the transfer line. The valve opening 
time-series are embedded into the code in the form of tables. The opening 
at arbitrary time is interpolated using these tables. 

The comparison of the actual opening with the interpolated values 
is shown in the Fig. 13. It can be seen from the figure that the dump 
valves and the line valves are operating according to their actual nominal 
behavior. The effect of the dump valve opening and closing on the liquid 
flow was verified earlier (see Sec. 6.5 in [LuchDG-II] . We now describe 
the comparison between the time-series data and model predictions for 
the pressure and temperature. 
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Figure 13. Relative valves opening for the valves cvll2, cvll7, cvllS, 
cvl36. The experimental time-series (black line) as compared to the 
interpolated values (red dashed lines). 


8.3 Pressure 

During the simulations that are discussed below the following simplifi- 
cations were introduced in the correlation module. No heat exchange 
was allowed between the phases. The heat exchange was between the 
gas and the wall, the liquid and the wall, and the environment and the 
wall was included. The rates of vaporization and the heat exchange are 
taken as constants. The flow is modeled as conceptually stratified (i.e. 
the flow regime is assumed to be stratified flow for all parameters). 

The simulations of the first 1600 sec of the chilldown process is shown 
in the Fig. 14 in comparison with the experimental pressure time-series 
data. It can be seen from the figure that the model can reproduce quite 
accurately the pressure drop along the pipe and the pressure oscillations 
during opening of the input valve MV151 and the main transfer valve 
Roll5. It can also reproduce the pressure jump at 195 sec at the location 
of the sensor PT161 and the distribution of the pressure alone the system 
(cf the change of the absolute pressure alone the transfer line at the 
locations PT PT102, PT157, PT161, PT148 in the figure). 

The simplified character of the correlations can be inferred by notic- 
ing the overestimation of the pressure level at the location PT161 by 16% 
for t > 500 sec and significant overestimation of the pressure oscillation 
at location PT148 at around 500 sec. These deviations are attributed to 
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Figure 14. Experimental time-series of pressure (black lines) in com- 
parison with model predictions (red dashed lines) for 6 sensors: PT102, 
PT157, PT161, PT148. 


the two main factors. 

Firstly, the only two-phase flow regime recognized in the simplified 
version is stratified flow. As a result, the liquid propagates quickly to the 
end of the transfer line while being in the contact with the hot wall. This 
leads to overestimation of the evaporation rates in the control volumes 
closer to the end of the line. Secondly, the heat transfer rate in the simpli- 
fied version is taken to be a constant approximately equal to the average 
heat transfer rate in the film boiling regime. This leads to an additional 
overestimation of the heat transfer rate. The overestimation is especially 
strong closer to the end of the transfer line where dryout, dispersed, and 
mist flow regimes are expected to dominate during chilldown. The over- 
estimation of the heat transfer leads to the overestimation of the mass 
flow and as a result to an increased pressure level in the end of the line. 
This simplifications will be lifted in the next version of the algorithm as 
will be described in the next report [LuchDG-III]. 

8.4 Temperature 

The model predictions for the temperature are compared with the exper- 
imental time-series data in the Fig. 15. A good agreement between the 
simulated and measured data can be seen un the figure. In particular, 
a good agreement observed at the location of sensor TT202 during first 
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200 sec is a substantial improvement as compared to the earlier results 
(see Sec. of [LuchDG-II]). 

This improvement was achieved by taking into account an additional 
section of the pipe that is located before valve MV151 and exposed to 
direct radiation. The corresponding changes included into the model 
geometry allowed for a better agreement with the experimental data not 
only at the location of sensor TT202, but also at all other locations as 
can be seen in the figure. 

However, the simplified character of the correlations used at this stage 
of the research do not allow to model a number of important features of 
the chilldown process. For example, strong oscillations of the temper- 
ature observed during cryogenic loading at KSC testbed could not be 
accurately reproduced in simulations. 

There are two main problems already mentioned at the and of the 
previous section. The hrst one is related to the fact that the current 
version recognizes only three flow regimes: (i) pure liquid; (ii) pure gas; 
and (iii) stratified flow. Second problem is that heat transfer coefficients 
are considered to be constants in the whole temperature range of the gas 
and liquid flow. 

As a result the strong peak in the heat transfer corresponding to 
the nucleate and transition boiling regimes (see Fig. 9) is not taken into 




t t 




t t 


Figure 15. Experimental time-series of temperature (black lines) in com- 
parison with model predictions (red dashed lines) for 6 sensors: TT202, 
TT105, TT162, TT165. 
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account. But it is this peak in heat (and accordingly in the mass) transfer 
rates that is responsible for the fast temperature drop observed in the 
experiment at approximately 167 sec in the top two figures and at 190 
and 395 sec in the bottom left figure of the Fig. 15. Such a fast cooling 
due to extensive evaporation can also cause a significant vapor lock with 
subsequent temperature raise as can be seen in the figure. To reproduce 
these features in the simulations a more accurate correlations that take 
into account the details of the boiling cure are required. 

Another deviation of the model predictions from the experimental 
data can be observed at the location of the sensor TT165. It can be seen 
from the figure that the model predicts faster temperature drop than 
observed during the test at around 500 sec. This result indicates the fact 
that the simplified version of the correlation module can only recognized 
stratified flow. As a result the heat transfer rate is overestimated in 
dispersed flow regimes observed experimentally at this location causing 
a faster temperature drop as compared to the experimental results (see 
temperature time-traces between 420 and 1000 sec in the bottom right 
figure of the Fig. 9). To avoid this problem the dispersed flow regime 
will be included in the next version of the code. 

Further details of the correlation module and the effect of the corre- 
lation parameters on the accuracy of the model predictions will be given 
in a separate report [LuchDG-III]. 

9 Conclusions 

In this report we discussed the details of the nearly-implicit method 
as it was codded in the C-|— |- version of the algorithm. We described 
calculations performed in the following main modules of the algorithm: 
(i) geometry; (ii) boundary conditions; (ii) correlations; (hi) first step, 
including solution of the expanded equations and of the velocity matrix; 
(iv) second step; and (v) control module. 

The following new features were added to the code as compared to 
the version discussed in the first two reports [LuchDG-I,LuchDG-II]: (i) 
the ability to integrate two-phase flow on an arbitrary non-uniform ID 
grid; (ii) improved module structure of the code that allows for more 
efficient computation and extends flexibility if the code with respect to 
future modifications; (iii) reduced integration time was achieved due to 
translation to C-|— extensive use of tables and linear interpolations, 
explicit codding of the matrix inverse for expanded conservation equa- 
tion and improved efficiency of the memory management; currently the 
integration is up to 400 times faster than real time; (iv) updated the cor- 
relation module. As a result the code is faster, more stable and allows 
for more accurate model predictions of the pressure and temperature 
dynamics. 

Substantial acceleration and improved accuracy of the C-|— |- version 
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of the algorithm paves the way for the future development of the on-line 
fault detection, isolation and recovery methods of loading control. These 
features make this method attractive for applications to autonomous 
control of cryogenic loading. 

However, the comparison of the results of simulations with the ex- 
perimental data also reveals some limitations of the current version. In 
particular, it was shown that temperature oscillations and vapor lock 
could not be accurately reproduced at some locations during initial stage 
of chilldown. Our analysis has shown that these limitations can be al- 
leviated by extending the correlation module. It was suggested that 
the next version of the correlation module should recognize more flow 
regimes of the two-phase flow and model more accurately the boiling 
curve under various flow conditions. The results of these improvements 
will be discussed in a separate technical report [LuchDG-III]. 


Nomenclature 

Acronyms 

Ax length of control volume 

q heat flux 

A pipe cross-section area 

b[ height of the liquid level in the pipe 
D pipe diameter 

E total specific energy 

e internal specific energy 

F friction factor 

/ friction coefficient 

H specific enthalpy 

h heat transfer coefficient 

hc,hfc convection and forced convection heat transfer coefficients 
Hig latent heat transfer of evaporation 
I perimeter of control volume 

My coefficient of the virtual mass 

p pressure 

Pr Prandtl number 
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Re Reynolds number 

S wall area of control volume 

T temperature 

u velocity 

2 : height of control volume center above the ground 

Greek Symbols 

a gas void fraction 

/3 liquid void fraction 

Xtt Lockhart-Martinelli parameter 

r mass flux 

K thermal conductivity 

fi viscosity 

p density 

a surface tension 

T drag force 

Superscripts 

^ upwind value 

'll) face centered value 

n index for the time step 

Subscripts 

, t time partial derivative 

, X spatial partial derivative 

dv dump valve 

9 gas 

j index for the staggered grid 

L index for the main grid 

I liquid 

s saturation value 
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